Caso Práctico 05: Colgate vs Crest

Master en Data Science para Finanzas: Predicción

library(openxlsx)
data <- read.xlsx("data.xlsx", colNames = TRUE, detectDates = TRUE)
head(data)

tail(data)

colgate <- ts(data$Colgate, start = 1958, frequency = 52)
crest <- ts(data$Crest, start = 1958, frequency = 52)

mydata <- ts(data[,c(3,4)], start = 1958, frequency = 52)

library(forecast)
library(ggplot2)
library(dplyr)

autoplot(colgate) +
  ggtitle("Colgate Market Shares") +
  xlab("Week") +
  ylab("Market Share")

autoplot(crest) +
  ggtitle("Crest Market Shares") +
  xlab("Week") +
  ylab("Market Share")

autoplot(mydata, facets = TRUE)

autoplot(mydata, facets = FALSE) +
  ggtitle("Cuotas de mercado de Crest y Colgate") +
  xlab("Tiempo") +
  ylab("Cuota de mercado")

ggseasonplot(colgate, month.labels = TRUE, month.labels.left = TRUE) +
  ylab("Market Share") +
  ggtitle("Seasonal plot: Colgate Market Share")

ggseasonplot(crest, week.labels = TRUE, week.labels.left = TRUE) +
  ylab("Market Share") +
  ggtitle("Seasonal plot: Crest Market Share")

ggseasonplot(colgate, polar = TRUE) +
  ylab("Market Share") +
  ggtitle("Seasonal plot: Colgate Market Share")

ggseasonplot(crest, polar = TRUE) +
  ylab("Market Share") +
  ggtitle("Seasonal plot: Crest Market Share")

ggsubseriesplot(colgate) +
  ylab("$ million") +
  ggtitle("Seasonal subseries plot: Colgate Market Share")

ggsubseriesplot(crest) +
  ylab("$ million") +
  ggtitle("Seasonal subseries plot: Crest Market Share")

qplot(crest, colgate) +
  ylab("Colgate market share)") + xlab("Crest market share")

gglagplot(crest)

gglagplot(colgate)

ndiffs(colgate)

[1] 1

ndiffs(crest)

[1] 1

nsdiffs(colgate)

[1] 0

nsdiffs(crest)

[1] 0

Fitting the ARIMA models

(fit_colgate <- auto.arima(colgate, seasonal = TRUE))

Series: colgate 
ARIMA(0,1,1)(1,0,0)[52] 

Coefficients:
          ma1    sar1
      -0.7618  0.0312
s.e.   0.0450  0.0658

sigma^2 estimated as 0.002281:  log likelihood=446.79
AIC=-887.59   AICc=-887.5   BIC=-876.74

(fit_crest <- auto.arima(crest, seasonal = TRUE))

Series: crest 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.6641
s.e.   0.0431

sigma^2 estimated as 0.002082:  log likelihood=458.98
AIC=-913.97   AICc=-913.92   BIC=-906.73

checkresiduals(fit_colgate)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)(1,0,0)[52]
Q* = 47.41, df = 53, p-value = 0.6907

Model df: 2.   Total lags used: 55

checkresiduals(fit_crest)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)
Q* = 40.09, df = 54, p-value = 0.9206

Model df: 1.   Total lags used: 55

fit_colgate %>% forecast(h = 16) %>% autoplot()

fit_crest %>% forecast(h = 16) %>% autoplot()

library(tsoutliers)
(colgate_outlier <- tso(colgate, types = c("TC", "AO", "LS", "IO", "SLS")))

Series: colgate 
Regression with ARIMA(0,1,1)(1,0,0)[52] errors 

Coefficients:
          ma1    sar1     TC43    AO102    LS136   TC196
      -0.8593  0.0084  -0.1298  -0.1606  -0.0997  0.1412
s.e.   0.0323  0.0674   0.0345   0.0416   0.0223  0.0344

sigma^2 estimated as 0.001902:  log likelihood=473.57
AIC=-933.15   AICc=-932.73   BIC=-907.83

Outliers:

plot(colgate_outlier)

(crest_outlier <- tso(crest, types = c("TC", "AO", "LS", "IO", "SLS")))

Series: crest 
Regression with ARIMA(0,1,1) errors 

Coefficients:
          ma1   LS136    AO167    TC196
      -0.7885  0.1646  -0.1427  -0.1322
s.e.   0.0447  0.0263   0.0392   0.0343

sigma^2 estimated as 0.001747:  log likelihood=484.43
AIC=-958.85   AICc=-958.63   BIC=-940.77

Outliers:

plot(crest_outlier)

summary(crest_outlier)

         Length Class          Mode   
outliers   5    data.frame     list   
y        276    ts             numeric
yadj     276    ts             numeric
cval       1    -none-         numeric
fit       19    forecast_ARIMA list   
effects  276    ts             numeric
times      3    -none-         numeric

(outliers_colgate_idx <- colgate_outlier$outliers$ind)

[1]  43 102 136 196

(outliers_crest_idx <- crest_outlier$outliers$ind)

[1] 136 167 196

(fit_colgate_outlier <- auto.arima(colgate_outlier$yadj, seasonal = TRUE))

Series: colgate_outlier$yadj 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.8598
s.e.   0.0312

sigma^2 estimated as 0.001867:  log likelihood=473.57
AIC=-943.13   AICc=-943.09   BIC=-935.9

(fit_crest_outlier <- auto.arima(crest_outlier$yadj, seasonal = TRUE))

Series: crest_outlier$yadj 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.7897
s.e.   0.0430

sigma^2 estimated as 0.001728:  log likelihood=484.43
AIC=-964.86   AICc=-964.81   BIC=-957.62

checkresiduals(fit_colgate_outlier)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)
Q* = 44.403, df = 54, p-value = 0.8211

Model df: 1.   Total lags used: 55

checkresiduals(fit_crest_outlier)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)
Q* = 35.724, df = 54, p-value = 0.974

Model df: 1.   Total lags used: 55

fit_colgate_outlier %>% forecast(h = 16) %>% autoplot()

fit_crest_outlier %>% forecast(h = 16) %>% autoplot()

Intervention Model

library(TSA)

arimax_model_crest <- arimax(crest,
                       order = c(0,1,1), #orden del ARIMA
                       xtransf = data.frame(I1 = (1*(seq(crest) >= 136))),             
                       transfer = list(c(1,0)), 
                       method = 'ML')

possible convergence problem: optim gave code=1

summary(arimax_model_crest)

Call:
arimax(x = crest, order = c(0, 1, 1), method = "ML", xtransf = data.frame(I1 = (1 * 
    (seq(crest) >= 136))), transfer = list(c(1, 0)))

Coefficients:

Se han producido NaNs

          ma1  I1-AR1  I1-MA0
      -0.7775  0.0401  0.1591
s.e.   0.0437     NaN     NaN

sigma^2 estimated as 0.001896:  log likelihood = 471.15,  aic = -936.31

Training set error measures:
                     ME       RMSE        MAE       MPE     MAPE
Training set 0.00205205 0.04346805 0.03385271 -3.397657 16.76844
                  MASE        ACF1
Training set 0.7955474 -0.01422081

autoplot(fitted(arimax_model_crest))+autolayer(crest)

checkresiduals(arimax_model_crest)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)
Q* = 45.428, df = 52, p-value = 0.7283

Model df: 3.   Total lags used: 55

escalon=(I1 = (1*(seq(crest) >= 136)))
w0=as.numeric(arimax_model_crest$coef["I1-MA0"])
d1=as.numeric(arimax_model_crest$coef["I1-AR1"])


plot(ts(stats::filter(escalon,filter=d1,method='recursive',side=1)*w0,
frequency=52,start=1958),type='h',ylab='')
abline(h=0)

arimax_model_colgate <- arimax(colgate,
                       order = c(0,1,1), #orden del ARIMA
                       xtransf = data.frame(I1 = (1*(seq(colgate) >= 136))),             
                       transfer = list(c(1,0)), 
                       method = 'ML')

summary(arimax_model_colgate)

Call:
arimax(x = colgate, order = c(0, 1, 1), method = "ML", xtransf = data.frame(I1 = (1 * 
    (seq(colgate) >= 136))), transfer = list(c(1, 0)))

Coefficients:
          ma1  I1-AR1   I1-MA0
      -0.8091  0.0006  -0.1031
s.e.   0.0412  0.5132   0.0495

sigma^2 estimated as 0.002163:  log likelihood = 453.02,  aic = -900.04

Training set error measures:
                       ME       RMSE        MAE       MPE     MAPE
Training set -0.001964486 0.04641939 0.03643914 -2.941924 12.55462
                  MASE       ACF1
Training set 0.7890364 0.05175557

autoplot(fitted(arimax_model_colgate))+autolayer(colgate)

checkresiduals(arimax_model_colgate)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)
Q* = 39.99, df = 52, p-value = 0.888

Model df: 3.   Total lags used: 55

escalon=(I1 = (1*(seq(colgate) >= 136)))
w0=as.numeric(arimax_model_colgate$coef["I1-MA0"])
d1=as.numeric(arimax_model_colgate$coef["I1-AR1"])


plot(ts(stats::filter(escalon,filter=d1,method='recursive',side=1)*w0,
frequency=52,start=1958),type='h',ylab='')
abline(h=0)

Funcion TRansferencia

mod.transf.colgate <- arimax(x = colgate, order = c(1, 0, 0), 
                                 xtransf = data.frame(crest), 
                                 transfer = list(c(0,20)), 
                                 method = "ML")
mod.transf.colgate

Call:
arimax(x = colgate, order = c(1, 0, 0), method = "ML", xtransf = data.frame(crest), 
    transfer = list(c(0, 20)))

Coefficients:
         ar1  intercept  crest-MA0  crest-MA1  crest-MA2  crest-MA3
      0.2349     0.4049    -0.4667    -0.0304    -0.0253     0.0337
s.e.  0.0618     0.0079     0.0546     0.0547     0.0565     0.0581
      crest-MA4  crest-MA5  crest-MA6  crest-MA7  crest-MA8
         0.0715    -0.0238     0.0695    -0.0197    -0.0058
s.e.     0.0583     0.0589     0.0591     0.0592     0.0595
      crest-MA9  crest-MA10  crest-MA11  crest-MA12  crest-MA13
        -0.0081     -0.0163      0.0121      0.0005      0.0619
s.e.     0.0596      0.0596      0.0606      0.0605      0.0607
      crest-MA14  crest-MA15  crest-MA16  crest-MA17  crest-MA18
          0.0214     -0.0533      0.0200      0.0308     -0.0517
s.e.      0.0607      0.0605      0.0599      0.0594      0.0583
      crest-MA19  crest-MA20
          -0.090      0.0954
s.e.       0.057      0.0567

sigma^2 estimated as 0.00157:  log likelihood = 463.15,  aic = -880.31

mod.transf.colgate <- arimax(x = colgate, order = c(1, 0, 0), 
                                 xtransf = data.frame(crest), 
                                 transfer = list(c(0,0)), 
                                 method = "ML")
mod.transf.colgate

Call:
arimax(x = colgate, order = c(1, 0, 0), method = "ML", xtransf = data.frame(crest), 
    transfer = list(c(0, 0)))

Coefficients:
         ar1  intercept  crest-MA0
      0.2870     0.4198    -0.4203
s.e.  0.0582     0.0077     0.0268

sigma^2 estimated as 0.0018:  log likelihood = 480.47,  aic = -954.94

library(lmtest)
coeftest(mod.transf.colgate)

z test of coefficients:

            Estimate Std. Error  z value  Pr(>|z|)    
ar1        0.2870442  0.0581593   4.9355 7.995e-07 ***
intercept  0.4198273  0.0077222  54.3665 < 2.2e-16 ***
crest-MA0 -0.4202765  0.0268355 -15.6612 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(mod1 <- Arima(colgate, order = c(1, 0, 0), xreg = crest))

Series: colgate 
Regression with ARIMA(1,0,0) errors 

Coefficients:
         ar1  intercept     xreg
      0.2870     0.4198  -0.4204
s.e.  0.0581     0.0077   0.0268

sigma^2 estimated as 0.00182:  log likelihood=480.47
AIC=-952.94   AICc=-952.79   BIC=-938.46

xnew <- fit_crest %>% forecast(h = 8)
mod1 %>% forecast(h = 8,xreg=xnew$mean) %>% autoplot()